Wannier-based calculation of the orbital magnetization in crystals 
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We present a first-principles scheme that allows the orbital magnetization of a magnetic crystal to 
be evaluated accurately and efficiently even in the presence of complex Fermi surfaces. Starting from 
an initial electronic-structure calculation with a coarse ab initio k-point mesh, maximally localized 
Wannier functions are constructed and used to interpolate the necessary k-space quantities on a fine 
mesh, in parallel to a previously-developed formalism for the anomalous Hall conductivity [X. Wang, 
J. Yates, I. Souza, and D. Vanderbilt, Phys. Rev. B 74, 195118 (2006)]. We formulate our new 
approach in a manifestly gauge- invariant manner, expressing the orbital magnetization in terms of 
traces over matrices in Wannier space. Since only a few (e.g., of the order of 20) Wannier functions 
are typically needed to describe the occupied and partially occupied bands, these Wannier matrices 
are small, which makes the interpolation itself very efficient. The method has been used to calculate 
the orbital magnetization of bcc Fe, hep Co, and fee Ni. Unlike an approximate calculation based on 
integrating orbital currents inside atomic spheres, our results nicely reproduce the experimentally 
measured ordering of the orbital magnetization in these three materials. 

PACS numbers: 71.15.-m, 71.15.Dx, 75.50.Bb 
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I. INTRODUCTION 

Magnetism in matter originates from two distinct 
sources, namely, the spin and the orbital degrees of free- 
dom of the electrons. In many bulk ferromagnets the 
spin contribution dominates, and it is therefore not sur- 
prising that the description of spin magnetism using first- 
principles methods is considerably more developed than 
that of orbital magnetism. In particular, the local spin- 
density approximation has been successful in studying 
magnetic materials for decades.^ 

Although the orbital moments in bulk solids are 
strongly quenched by the crystal field and typically give 
small contributions to the net magnetization — between 
5% and 10% in Fe, Co, and Ni — they can be mea- 
sured very accurately with the help of gyromagnetic 
experiments.^''^ Moreover, there are known instances 
where orbital magnetism plays a prominent role. Some 
examples include weak ferromagnets with large but op- 
posing spin and orbital moments, low-coordination 
systems such as magnetic nanowires,^ and the recently- 
predicted large orbital magnetoelectric coupling in topo- 
logical insulators and related materials. In addition, 
magnetic resonance parameters such as the NMR^^"^"^ 
and EPR^** shielding tensors can be conveniently cal- 
culated as the change in (orbital) magnetization under 
appropriate perturbations. These examples highlight 
the need to develop accurate and efficient first-principles 
schemes for describing orbital magnetism in solids. 

The traditional way of computing the orbital mag- 
netization M is by integrating currents inside atom- 
centered muffin-tin spheres. ^^'^^ This requires choosing. 



somewhat arbitrarily, a cutoff radius and neglects con- 
tributions from the interstitial regions. A rigorous the- 
ory for the orbital magnetization of periodic crystals 
free from such uncontrolled approximations was obtained 
only recently^^"^° (see Ref. 21 for a review). The theory 
was developed in the independent-particle framework, 
and its central result is the expression 



M = ^ E / 1'^^] 



(1) 



(9kUnk| X (^^k + £nk - 2£_F)|c?kWnk) 



Here the integral is over the Brillouin zone (BZ), [dk] 
stands for dk/{2Trf, i7k = e-*-''?ie*-'" is the Bloch 
Hamiltonian whose eigenfunctions are the cell-periodic 
Bloch functions |unk) with eigenvalues e„k, /nk is the 
zero-temperature Fermi occupation factor, and ep is the 
Fermi energy. The third term in Eq. (1) vanishes in or- 
dinary insulators, but must be included in the case of 
metals. ^^"2° 

The implementation of Eq. (1) requires a knowledge 
of the k-space gradients |9kW„k) of the occupied Bloch 
states. An easier quantity to compute in practice is 
the covariant derivative |9k'"nk), defined as the projec- 
tion of |(?kU„k) onto the unoccupied bands. It turns out 
that the replacement \du_Unk) — >■ |9kUnk) in Eq. (1) leaves 
the sum of its terms invariant. For band insulators the 
covariant derivative can be conveniently evaluated by fi- 
nite differences.^® Unfortunately, the discretized covari- 
ant derivative approach cannot be applied to metals, as 
it relies on having a constant number of occupied bands 
throughout the BZ. 
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Thus far, the only first-principles applieation of Eq. (1) 
to metals is the calculation in Ref. 14 of the spontaneous 
orbital magnetization in Fe, Co, and Ni crystals, where 
the k-derivatives of the Bloch wave functions were evalu- 
ated using a linear-response method. This carries a cost 
per k-point comparable to that of a non-self-consistent 
ground-state calculation. The number of k-points needed 
to converge the BZ integral in Eq. (1) for Fe, Co, and Ni 
is quite significant, rendering the full calculation rather 
time-consuming. 

In this work we develop an alternative approach which 
greatly reduces the computational cost of evaluating 
Eq. (1) for metals. Our implementation relies on a 
method for constructing well-localized crystalline Wan- 
nier functions (WFs) by post-processing a conventional 
band structure calculation. ■^■^ A key ingredient is the 
"band disentanglement" procedure of Ref. 24, which al- 
lows one to obtain a set of WFs spanning a space that 
contains the occupied valence bands as a subspace. These 
WFs are essentially an exact tight-binding basis for those 
ab initio bands that carry the information about M. 
Working in the Wannier representation, the problem of 
evaluating M can then be reformulated in a very eco- 
nomical way. This reformulation involves setting up the 
matrix elements of certain operators in the Wannier ba- 
sis. Once that is done, the integrand of Eq. (1) can be 
evaluated very inexpensively and accurately at arbitrary 
points in the BZ. The cost per k-point of the entire pro- 
cedure is significantly reduced, especially in cases where 
a dense sampling of the BZ is needed to achieve conver- 
gence. Our method builds on the work of Ref. 25, where a 
similar "Wannier interpolation" strategy was introduced 
to calculate the anomalous Hall conductivity (AHC) of 
ferromagnetic metals. 

This manuscript is organized as follows. In Sec. II 
the orbital magnetization formula, Eq. (1), is recast in 
a gauge-invariant form, and a related expression for the 
AHC is introduced. We then describe step by step the 
formalism used to express physical quantities in the Wan- 
nier representation. That formalism is applied in Sec. Ill 
first to the AHC, and then to the orbital magnetization. 
In Sec. IV we describe the procedure for evaluating the re- 
quired k-space matrices by Fourier interpolation. Some 
details of the first-principles calculation and Wannier- 
function construction are given in Sec. V, followed by an 
application of the method to bcc Fe, hep Co, and fee Ni 
in Sec. VI. We conclude in Sec. VII with a brief summary 
and discussion. 



II. PRELIMINARIES 



Orbital magnetization and anomalous Hall 
conductivity 



For our purposes it will be convenient to recast Eq. (1) 
in a different form as introduced in Ref. 19. We begin by 



writing the axial vector M as an antisymmetric tensor, 



Mal3 = -Cafi'^M^ , 



(2) 



where Greek indices denote Cartesian directions. We now 
partition M^p into two terms. 



The "local circulation" is 



(3) 



and the "itinerant circulation" is 

The k-dependent quantities Ga/3, and Ha/3 are 



= Tr [{daP)Q{dfsP) 
Go^p = Tr [{daP)QHQ{dfiP) 
H^^ = Tr \H{daP)Q{d0P) 



(6) 
(7) 
(8) 



where "Tr" denotes the electronic trace, 9„ stands for 
d/dka, and the subscript k is implied in the Bloch 
Hamiltonian H and in the projection operators P and 
Q = 1—P spanning the occupied and unoccupied spaces 
respectively (here 1 is the identity operator in the full 
Hilbert space). 

We work at T = 0, so that Eqs. (2)-(8) yield the same 
result as Eq. (1). This can be seen by writing P in terms 
of the Bloch eigenstatcs. 



(9) 



and setting the occupancies /„ to either one or zero. 

Compared to Eq. (1), the above formulation has 
the advantage of being manifestly gauge-invariant, i.e., 
independent of any k-dependent phase twists applied 
to the occupied Bloch states, or more generally, any 
k-dependent unitary mixing among them. Because 
Eqs. (6) (8) are written as traces involving projection 
operators, they remain valid no matter how we choose to 
represent the occupied space at each k. Instead, Eq. (1) 
is written explicitly in terms of the energy eigenstates, 
that is, it assumes a Hamiltonian gauge. 

It should come as no surprise that it is possible to 
cast a physical observable such as the orbital magne- 
tization in a gauge-invariant form. More interestingly, 
the two terms in Eq. (3) are individually gauge-invariant, 
and this led to the speculation that they might be sepa- 
rately observable. That is indeed the case — at least in 
principle — as discussed in Ref. 26. 
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Before continuing wc mention another physical observ- 
able, the intrinsic anomalous Hall conductivity, which 
can be expressed in gauge-invariant form as 

''^0=-J jmi-21mF„0) . (10) 

In the Hamiltonian gauge the integrand of this equation 
acquires a more familiar form, i.e., as the k-space Berry 
curvature summed over the occupied bands. "'^^ 

While developing the formalism in Sec. Ill we shall 
first consider the AHC before tackling the more complex 
case of the orbital magnetization. This will allow us to 
make contact with the work of Ref. 25, where a Wannier 
interpolation scheme was developed for the AHC, but 
using a somewhat different formulation. 

B. Wannier space and gauge freedom 

The crux of our approach is to express the observables 
of interest (orbital magnetization and AHC) not in terms 
of the Bloch cigcnstatcs |urik), but using an alternative 
set of Bloch-like states constructed at each k as ap- 
propriately chosen linear combinations of energy eigen- 
states. The defining feature of the new states is that they 
are smooth functions of k. As a result, the corresponding 
Wannier functions (WFs) 

|I^") = 7^Ee"'''^''"'^l"nk) (11) 
k 

(here is the number of k-points distributed on a uni- 
form BZ mesh) are well-localized in real space, and for 
this reason we shall say that the states belong to the 
Wannier gauge. The ability to construct a short-ranged 
representation of the electronic structure in real space 
is what will allow us to devise an efficient and accurate 
interpolation scheme for the k-space quantities ImF„^, 
Im Gafs , and Im Hafs ■ 

The construction of a Wannier basis proceeds in two 
steps, which we call "space selection" and "gauge selec- 
tion." In the case of insulators, the space selection is 
typically obvious; we want the WFs to span just the oc- 
cupied subspacc. In k-space, we represent this subspace 
by a k-depcndent "band projection operator" Pk defined 
as in Eq. (9) (with /„ = 1 and for occupied and empty 
bands respectively). Denoting the Wannicr-space projec- 
tion operator by Pk, we can then set Pk = -Pk- 

For metals, on the other hand, the space selection step, 
also known as "band disentanglement," is more subtle. 
One wants to choose a J-dimensioiial manifold, repre- 
sented by the projection operator Pk, throughout the 
BZ such that it has the following properties: (i) it must 
contain as a subspacc the set of all occupied eigenstates 
(hence J cannot be less than the highest number of oc- 
cupied bands at any k); (ii) it must display a smooth 
variation with k, in the sense that Pk is a difi'erentiable 



function of k. A procedure for extracting an optimally- 
smooth space from a larger set of band states was devel- 
oped in Ref. 24. The resulting space typically contains 
some admixture of low-lying empty states in addition to 
the occupied states. 

The gauge selection consists of representing the 
smoothly-varying space Pk using a set of J Bloch-like 
states which are themselves smooth functions of k, 

I'k = El^^k)(Ckl- (12) 

n 

From these the WFs are constructed via Eq. (11). 

A procedure for selecting an optimally-smooth Wannier 
gauge was developed in Ref. 23, such that the result- 
ing WFs are maximally-localized in the sense of having 
the smallest possible quadratic spread. The localization 
procedure was originally devised with an isolated group 
of bands in mind (e.g., the valence bands of an insula- 
tor), but it can be applied to any smoothly- varying Bloch 
manifold of fixed dimension J. 

In the Wannier gauge the projected Hamiltonian Hk = 
IPk-f^klPk takes the form of a non-diagonal J x J matrix, 

H^„(k) = (n^kl^kl^rk) • (13) 

We define the Hamiltonian gauge in the projected space 
as the gauge in which this matrix becomes diagonal, 

<„(k) = C/t(k)HW(k)[/(k) = e„k5nm • (14) 

Because of the nature of the space selection step, the 
"projected eigenvalues" £„k agree with the true ab ini- 
tio eigenvalues e„k for all occupied states, but they may 
differ for unoccupied states. 

The unitary matrices U{\i) that diagonalize H^(k) can 
be used to transform other objects between the Wannier 
and Hamiltonian gauges. For example, the Bloch states 
transform as 

J 

hn^.) = E A)Umn{^) • (15) 

m 

The gauge-invariance of the projection operators Pk can 
be checked explicitly by inserting the right-hand-side of 
Eq. (15) in place of |u^) in Eq. (12). 

C. Projection operators and occupation matrices 

In Eq. (9) we introduced the projection operator Pk 
onto the occupied manifold at k, and in Eq. (12) the 
projection operator Pk onto the Wannier space at k. Fig- 
ure 1 represents schematically the relationship between 
those two subspaces, as well as other related subspaces 
to be defined shortly. The notation in the figure is as fol- 
lows: a double staff is used for objects that concern the 
distinction between the space spanned by the WFs ("in- 
ner") and the corresponding orthogonal space ("outer"). 
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FIG. 1. The Hilbert space of Bloch functions at wavevector k 
can be decomposed either as 1 = P + Q (occupied and unoc- 
cupied spaces) or as 1 = P + Q ("inner" and "outer" spaces, 
where the "inner" space is spanned by Wannier functions). 



while a single staff will be used for objects that distin- 
guish between the occupied and unoccupied parts of the 
Wannier space. 

Consider, for example, the projection operator P onto 
the occupied bands, as defined in Eq. (9). Recall that 
k labels were suppressed in Sec II A; with k temporarily 
restored, we have 



n 
J 

— ^ l^mk) fmn,k 



(16) 



nkl 1 



where ^ = (Wmkl-Pklw^) is the (non-diagonal) occu- 
pation matrix in the Wannier gauge. In the following, we 
will use a strongly condensed notation, leaving band in- 
dices (and sums over them) implicit, omitting wavevector 
subscripts, and dropping the superscript "W" from ob- 
jects in the Wannier gauge. So, for example, we will write 
the Wannier-gauge Bloch states simply as \u) instead of 



w 



In this notation, Fk is expressed in the Wannier gauge 
as just 



P = \u)f{u\ . 



(17) 



Similarly, the projector onto the Wannier ("inner") space 
is 



We also define 



with 



P = \u){u\ . 

Q = i -p , 

Qin = \u)g{u\ , 



(18) 

(19) 
(20) 

(21) 
(22) 



where "1" denotes the J x J identity matrix. 

In practice the occupation matrix is first evaluated in 
the Hamiltonian gauge in which it is diagonal, and then 
rotated to the Wannier gauge with the help of Eq. (15), 

/ = C//"C/^ . (23) 

The matrices / and g are idempotent, and satisfy 

f9 = gf = 0, (24) 



as well as 



where 



[f,M] = [g,M]=0 



H = {u\H\u) 



(25) 



(26) 



is Eq. (13) using the concise notation. Equations (24) 
and (25) imply 



mg = o. 



(27) 



Note that /H = H/ = /H/. We shall also make frequent 
use of relations such as QinQ = QQm = 0, = Q, 

Qin — Qin'j CtC. 



D. Compendium of "geometric" matrix objects 

We list here for future reference a number of addi- 
tional J X J matrices that will be used to express Im Fa/3, 
ImGap, and ImHa^ in the Wannier representation: 

(28) 
(29) 

- A„A0 , (30) 

- dfsAa , (31) 
i[Aa,Ap] . (32) 

The Hermitian matrices A^, fla/3, and Da^S sse known 
as the Berry connection. Berry curvature, and gauge- 
covariant Berry curvature. They play a central role in 
the theory of geometric-phase effects in solids. ^'^ 

In addition to the above objects, which represent in- 
trinsic properties of the Bloch manifold, we shall make 
use of a number of similarly-defined quantities which also 
depend on the Hamiltonian (and therefore are not strictly 
speaking "geometric"): 



A„ 


= i{u\dau) , 






= {daU\dfsu) , 






= {dau\il\di3u) 








= a„A 




= ^F„^ - iWl^ 





= i{u\H\dau) , 
Ba = i{u\H^\dau) = Ba- HA„ 

= {daU\H\di3u) , 

= - A„P^ - bIA^ , 
= Aafs-i (a^Bp - ApBa 



h.c. 



(33) 
(34) 
(35) 

(36) 
(37) 

(38) 
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where "h.c." stands for Hcrmitian conjugate. Note that 
while k„jj and A^,^ are Hermitian, IBq, is not. 

All of the matrices listed above are written in the Wan- 
nier gauge, where they are smooth functions of k and can 
be evaluated efficiently using Fourier transforms, as will 
be described in Sec. IV. 



III. WANNIER-SPACE REPRESENTATION OF 
PHYSICAL QUANTITIES 

A. Anomalous Hall conductivity 

1. Derivation 

As a first application of our formalism, let us consider 
Eq. (10) for the AHC. The integrand is the Berry cur- 
vature — 2ImFct/3, and we wish to write it as a trace of 
products of matrices defined within the Wannicr space. 

Our starting point is Eq. (6) for Fq,^. In preparation 

for taking the trace therein, let us express {do:P)Q in the 
Wannier gauge. Differentiating Eq. (17) leads to 



With this relation in hand it becomes straightforward 
to evaluate Eq. (6), and we quickly arrive at 



daP = \daU)f{u\ + \u)f{daU\ + \u)fa{u\ 



(39) 



where = daf. Multiplying on the right with Q and 
using Eq. (20) yields two terms for {daP)Q- One is 

(9„P)Q = \u)f{dc,u\(Q , (40) 

where = was used, and the other is 

{daP)Qin = |u)/(+iA„)5(u| + \u)Ug{u\ , (41) 

where Eq. (24) was used. Now, from Eq. (23) we find^^ 

./« = * [./, Ja] , (42) 

where the Hermitian matrix ,/„ , like / itself, is first eval- 
uated in the Hamiltonian gauge, being defined as 

J« = lU^d^U , (43) 

and then rotated into the Wannier gauge, 

J„ = UJ^U^ . (44) 

Using Eq. (42) in Eq. (41) and defining 

Ao, = Aa + Ja, (45) 

we arrive at the compact form 

{daP)Qin = i\u)fA^g{u\ . (46) 
The desired expression is the sum of Eqs. (40) and (46), 

{daP)Q = \u)f{dau\<Q + i\u)fAag{u\ . (47) 



(48) 



where "tr" denotes the trace over J x J matrices, not 
to be confused with the trace "Tr" over the full Hilbert 
space introduced earlier. We are interested in the imagi- 
nary part, and using Eq. (32) we find 

-2ImF„^ = Retr[/D„^] - 2Imtr[fA„gA0] (49) 
= Retr[/Cl„^] + 2Imtr[/A„A^ - fA^gA^] . 

Expanding Aa and using the identity 

Imtr[/A„A^] = Imtr[/A„(/ + g)Ap] (50) 
= Imtr[/Aa5A^] , 

we end up with 

-2ImF„^ = Retr[/n„/3] (51) 
- 2Imtr[/Aa.gJ/3 + fJagAp + fJagJii] ■ 

This expression for the Berry curvature of the occupied 
states is our first important result. 

2. Discussion 

The ingredients that go into Eq. (51) are the matrices 
Aq,, Oa^, /, and Ja expressed in a smooth (but otherwise 
arbitrary) gauge. It should be noted that while A^ and 
flap are themselves smooth functions of k, this is not so 
for / and Jq,. Consider /, given by Eq. (23). In metals 
it is affected by the step-like discontinuity in at the 
Fermi surface. More generally it is also affected by the 
wrinkles in the rotation matrix U (recall that U relates 
via Eq. (15) the smooth Wannier-gauge Bloch states to 
the Hamiltonian-gauge states, which are non-analytic as 
a function of k at points of degeneracy). The situation is 
even more severe in the case of Ja- Because it contains 
the derivative daU of a non-smooth function, it has spikes 
in k-space. 

How does one reconcile the existence of irregular and 
spiky quantities inside Eq. (51) with the form of Eq. (6), 
which suggests that Fafi is a smooth function of k, ex- 
cept possibly when crossing the Fermi surface (when a 
state comes in or out of the occupied manifold, intro- 
ducing a discontinuity in Pk)? The answer is that while 
Ja itself has a very irregular behavior, combinations like 
fJag which actually appear in Eq. (51) do not, as will 
be discussed in Sec. IV A. 

Let us now make contact with the formulation of Wan- 
nier interpolation developed in Ref. 25. In that work 
the expression for the Berry curvature of the occupied 
states [Eq. (32) therein] was written in the Hamiltonian 
gauge, where the occupation matrices are diagonal. This 
required transforming A„ and flap from the Wannier 
gauge, where they are first constructed, to the Hamil- 
tonian gauge. Here we choose to keep everything in the 
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Wannicr gauge throughout. The advantage is that even 
though the matrices / and Jq, have to be constructed first 
in the Hamihonian gauge, it is straightforward to rotate 
them into the Wannier gauge where all other needed ob- 
jects are constructed. Instead, the reverse transforma- 
tion of those other objects can in certain cases become 
nontrivial. 

The two formulations are of course equivalent, and it 
is instructive to recover explicitly from Eq. (51) the cor- 
responding expression in Ref. 25. Consider for example 
the last term in Eq. (51). Taking the trace in the Hamil- 
tonian gauge we find 



tr[fJ^gj0] = fnJlnmi^ ' f^)jf,mn , (52) 

nm 

and thus 



-2Imtr[/Jc,yJ/3] 



(53) 



which agrees with the last term in Eq. (32) of Ref. 25 
(Z)i^^ therein corresponds in our notation to —iJ^). 

It is pleasing to see that Eq. (51), when converted 
to the Hamiltonian gauge, reduces to what was termed 
in Ref. 25 the "sum over occupied bands" expression, 
where individual terms have spiky features only when 
two bands, one occupied the other empty, almost touch 
at the Fermi level, as the factor fn^fm suppresses spikes 
associated with pairs of occupied states. As we shall see 
in Sec. IV A, this is a general feature of our formulation, 
which leads naturally to expressions where the spiky ob- 
ject Ja appears under the trace sandwiched between / 
and g. 



B. Orbital magnetization 

Let us now apply to Im Hafj and Im Gap the same 
strategy developed above for Im Fap, completing the list 
of quantities needed to evaluate the orbital magnetiza- 
tion. We remark that it would be possible to proceed 
along the lines of Ref. 25 in order to arrive at "sum 
over occupied bands" expressions for those quantities in 
the Hamiltonian gauge. However, we found such an ap- 
proach to be rather cumbersome, especially in the case 
of ImGafi- 



1. Derivation 



Inserting Eq. (47) into Eq. (8) leads to 



so that 

-2lmHo,p = Retr [/H/O^^j] - 2Imtr [fUfA^gA^] 
= Retr [/H/n„^] 

+ 2 Im tr [/H/A„ - fUfA^gAp] . (55) 
Using Eq. (45), this takes the desired form 

-21m H^p = Retr [/H/n„^] + 2Imtr [fUfAafAp] 
- 2 Imtr [/H/ {A^gJ^ + J^gA^ + J^gM 

(56) 

in terms of basic matrix objects with every sand- 
wiched between an / and a g (after taking account of 

the cyclic property of the trace) . 
Repeating for Ga^, we find 



and 



tr 



/ + A^gB^ + BigA^ + A„gMgAp^ 



-2ImGo 



2Imtr 



Retr 



(57) 

(58) 



(^fAc^gB^ -a-^/3^+ fA^gMgA^ 



Expanding and B^, this becomes 

-2ImGc./3 =Retr [fA^p] 



2 Imtr 



/ {AagBp - A^gBAp - A^Bfi -a^, 



f (A„HA;3 + A^gUgA^) 



(59) 



Next we expand A^ and gather terms in three groups, 
containing zero, one, and two occurrences of the matrices 
Jce and Jg,^^ 

-2ImGap = JO + Jl + J2 . (60) 
The J2 group contains only one term, 

J2 = -2Imtr[/J„5H5J;3] . (61) 
The Jl group is 
Jl = -2 Imtr [fJag (B^ - MAp + MgA/s) - a o /3] 



-2 Imtr [/ Ja5lB;3 - a -H- 



(62) 



where in the second equality we replaced one instance 
of g with 1 — / and then used Eq. (27). The JO group 
reads, after combining certain terms and canceling out 
some others, 

JO = Retr[/A„;3] (63) 
+ 2Imtr[/(A„/B/3-A;3/B„-A„H/A^)] . 

This can be simplified further with the help of the fol- 
lowing identity proven in the Appendix, 



(54) 



tr IfAafBp] = tr [/A^/HA^j] , 



(64) 
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which leads to 

JO = Retr [/A„^] - 2 Imtr [/H/A„/A^] . (65) 

Collecting terms, we find 

-2ImGa/3 = Rctr[/A„^i] - 2 Imtr [/H/A„/A^] (66) 
- 2 Imtr [/ {Ja9^p - a ^ + fJ^gUgJ^] . 



2. Final expressions 

The quantities ImFap, ImG„^, and Imi?„/3 enter 

the orbital magnetization expression in the combinations 
lm.{Gai3 —spFafj) and Im(ifQ,/9 — EF-Fa/?)- Using the con- 
densed notations 



X" = fXf , 
X' = gXg , 
X+ = gXf , 
X- = fXg , 



(67) 
(68) 
(69) 
(70) 



in Eqs. (51), (56), and (66), we obtain for the integrand 
ofMic, Eq. (5), 

+ Retr [(hO - sf) 0°^] + 2 Imtr [H^AO A^] (71) 



2 Imtr 



(H" - sf) (a- J+ + J-A+ + J- J+) 



and for the integrand of M^^, Eq. (4), 

-2Im(Gc</3 - SFFap) = 

+ Retr [AO^ - £^^0°^] - 2 Imtr [hOaO A^] (72) 
-2Imtr|[j-(B^ -efAJ) - a ^ p] 
+ J-(Hi-£^)J+}. 

Note that the second terms in these two equations are 
equal and opposite. For the special case of an insulator 
with / = 1 and .g = 0, only the first two terms are nonzero 
in each of Eqs. (71-72), and these expressions reduce to 
those derived in Ref. 19. 



IV. INTERPOLATION OF THE 
WANNIER-GAUGE MATRICES 

We calculate the orbital magnetization by averag- 
ing Eqs. (71) and (72) over a sufficiently dense grid of 
k-points in order to approximate the BZ integrals in 
Eqs. (4) and (5). At each k the matrices /, J^, H, Ika, 
IBa, Oq^, and h^p are needed, and they are calculated 
by Fourier interpolation as follows. 



A. Fourier transform expressions 



We start with the matrix H. Inverting Eq. (11), 



R 

and inserting into Eq. (26) yields 



,ik-R, 



0|-H|R) . 



(73) 



(74) 



R 



We emphasize that any desired wavevector k can be 
plugged into this expression, allowing one to smoothly 
interpolate the matrix H between the ah initio grid 
points used in Eq. (11) to construct the WFs. Diagonal- 
izing H [Eq. (14)] and using the resulting rotation matrix 
U and interpolated eigenvalues e„ in Eq. (23) yields the 
occupation matrix /. 

Next we consider the matrix J a- In practice it can be 
calculated by inserting U into Eq. (44) and then using^^ 



(75) 







if n 



m 



where Ha = 9alH is obtained by differentiating Eq. (74). 
In the vicinity of band degeneracies and weak avoided 
crossings the denominator in Eq. (75) becomes small, 
leading to strong peaks in Janm ^ function of k. If 
both bands n and m are occupied, such peaks must even- 
tually cancel out in the final expressions for the AHC and 
orbital magnetization, as discussed in Sec. Ill A 2. 

We can make such cancellations explicit from the out- 
set by noting that Ja only appears in the combinations 
JqT — fJag and J+ = gJaf- Taking the former, for 
example, we find using Eqs. (23) and (44) that 

J- = C/J«-C/t (76) 
where the matrix elements of = f^J^ (l — /^) are 



J"- 



if n occupied and m empty, 
otherwise. 

. ! . . . . 

Janm givcn by a similar expression, but with m oc- 
cupied and n empty. Unlike Eq. (75), the expressions 
for and J^"*" are well behaved in that they will only 
show peaks when the direct gap is small and mixing be- 
tween occupied and empty states is strong. By working 
directly with them, we avoid introducing any quantity 
that would react strongly to crossings among occupied 
states. 

While /, J~ and J+ are first calculated in the Hamil- 
tonian gauge and then converted to the Wannier gauge, 
the remaining quantities entering Eqs. (71) and (72) are 
most easily calculated directly in the Wannier gauge, in 
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the same way as H. It is sufficient to consider the Wan- 
nier representation of the three basic quantities Aq., IBq., 
and <Ca^ introduced in Sec. II D. Inserting Eq. (73) into 
the respective definitions we find 



A„=5^e''^-^(0|r„|R), (78) 

R 

B„=^e^''-^(0|?i(f-ii!)„|R) , (79) 

R 

<Ca0 = e'^''^{0\ran{f - -R)/3|R) • (80) 

R 

The expressions for flap and Aap are obtained by insert- 
ing Eq. (78) into Eq. (31) and Eq. (80) into Eq. (37), 
respectively. 

It was shown in Ref. 25 that the computation of the 
AHC requires a knowledge of the Wannier matrix ele- 
ments of H and r. Inspection of the Fourier transform 
expression given above reveals that the bulk orbital mag- 
netization requires in addition the matrix elements of "Hr 
and rUr. This is more than might have been antici- 
pated, given that the matrix elements of r and Hr are 
not needed for calculating the orbital moment of finite 
samples under open boundary conditions, but it is the 
price to be paid for a formulation that extends also to 
the case of periodic boundary conditions. 



B. Evaluation of the real-space matrices 

We shall follow the approach of Ref. 25, whereby the 
needed real-space matrix elements are actually evalu- 
ated in reciprocal space. Inverting the Fourier sums in 
Eqs. (74), (78), (79), and (80), we find 



(0|?i|R) = ^5^e-*-^Hk 



(81) 



(0|f«|R) = e-^'"^(wk|a««k) 



e-*-^w;b6„Mkb , (82) 



k,b 



{0\n{f - RUK) = ^J2 e-^'^'Kl^kiaaWk) 



_±_ e-*-^«;b6„Hk,b , (83) 



k,b 



- 7^ *'""'W'6i&laW'62&2^Hk,bi,b2 



k,bi,b2 



(84) 



where the sums are over the A/"^ points in the uniform 
ab initio mesh. The second equalities in Eqs. (82)-(84) 



follow from using a finite-difference expression for the 
derivatives of the smooth Bloch states, 



\daUk} ='Ywbba\uk+b) + 0{b'^) 



(85) 



{wb are appropriately chosen weights, and the sum is over 
shells of vectors b connecting a point k on the ab initio 
grid to its neighbors), together with the definitions 

Hk,b = (wkl-^kl^k-Kb) , (86) 
Mk,b = (wkluk-Hb) , (87) 
Hk,bi,b2 = {uk+hi\Hk\uk+h2) I (88) 

which complement Eq. (26) for Hk, needed in Eq. (81). 
Writing the states \u) as linear combinations of the orig- 
inal ab initio eigenstates lu*^). 



\Unk) = J2 l"mk>^k,, 



(89) 



we arrive at the following expressions for the J x J ma- 
trices appearing in Eqs. (81)- (84): 



Hk = FjHkVk , 
Hk,b - T/J?{kA4k.bM.+b , 
Mk,b - Vj>fkVk+b , 

Hk,bi,b2 = V^k+bi^k.bi.bsVk+ba 



where 



-Hw = «|i^k|< 
Xk.b = «K+b) 

Wk,bi,b2 



^k+b 



(90) 
(91) 
(92) 
(93) 

(94) 
(95) 
(96) 



The diagonal eigenvalue matrix "Hk and the overlap 
matrix A^k,b are readily available, as they constitute the 
input to the space-selection and gauge-selection steps in 
the wannierization procedure. While they suffice for cal- 
culating the AHC,^^ as well as M^^, the M'"^ term in the 
orbital magnetization requires the additional quantities 

^k,bi,b2 ■ 



V. COMPUTATIONAL DETAILS 

Planewave pseudopotential calculations were carried 
out for the ferromagnetic transition metals bcc Fe, 
hep Co, and fee Ni at their experimental lattice 
constants (5.42, 4.73, and 6.65 bohr, respectively). 
The calculations were performed in a noncollinear 
spin framework, using fully-rclativistic norm-conserving 
pseudopotentials"^° generated from similar parameters as 
in Ref. 25. The energy cutoff for the expansion of the 
valence wavefunctions was set at 120 Ry for Fe and Ni 
and at 140 Ry for Co; a cutoff of 800 Ry was used 
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for the charge density. Exchange and correlation ef- 
fects were treated within the PBE generahzed-gradient 
approximation . ^ ^ 

The calculation of the orbital magnetization comprises 
the following sequence of steps: (i) self-consistent total- 
energy calculation; (ii) non-self-consistent band structure 
calculation including several conduction bands; (iii) eval- 
uation of the matrix elements in Eqs. (95) and (96); (iv) 
wannierization of the selected bands; and (v) Wannier 
interpolation of Eqs. (71) and (72) across a dense k- 
point mesh, with the value of e_F taken from step (i). 
Steps (i) and (ii) were carried out using the PWSCF 
code from the Quantum-Espresso package,^^ and in 
step (iii) we used the interface routine pw2wannier90 
from the same package, modified to calculate Eq. (96) in 
addition to Eq. (95). Step (iv) was done using the Wan- 
NIER90 code,'^^ and for step (v) a new set of routines was 
written (we plan to incorporate these in a future release 
of the WANNIER90 distribution). 

The BZ integration in step (i) was carried out on 
a 16x16x16 Monkhorst-Pack mesh,^^ using a Fermi 
smearing of 0.02 Ry. In step (ii), the 28 lowest band 
states were calculated for bcc Fe and fee Ni on A/" x A/" x A/" 
k-point grids including the P-point. (For hep Co, with 
two atoms per cell, 48 states per k-point were calculated.) 
After testing several grid densities for convergence (see 
Sec. VIA below), we settled on A/" = 10 for all three ma- 
terials. In step (iv) we followed the procedure described 
in Ref. 35 to generate eighteen disentangled spinor WFs 
per atom, capturing the s, p, and d characters of the 
selected bands. In the case of fee Ni, we also tested 
an alternative set consisting of only fourteen WFs, ten 
of which are atom-centered d-like orbitals while the re- 
maining four are ,s-like and are centered at the tetra- 
hedral interstitial sites. We found excellent agreement 
- to within 0.0002 /zs/atom - between the values of M 
obtained with the two sets of WFs. 

It should be kept in mind that our calculations use 
a pscudopotential framework in which the contributions 
to M coming from the core region arc not quite de- 
scribed correctly. A rigorous treatment using the so- 
called GIPAW approach^^ was developed in Refs. 11 and 
14. It was shown that Eq. (1), written in terms of the 
pseudo-wavefunctions and pseudo-Hamiltonian, must be 
supplemented by certain core-reconstruction corrections 
(CRCs) in order to obtain the full orbital magnetization. 

We know from the work of Ref. 14 that the CRCs are 
small for bulk Fe, Co, and Ni, of the order of 5%. This 
suggests that the errors inherent in our uncontrolled ap- 
proximations in the core are also of the same order. If 
one wants to treat the problem correctly and capture this 
missing 5%, one should use the GIPAW approach. How- 
ever, the issues of implementing the CRCs are completely 
orthogonal to the issues of Wannier interpolation, and so 
we have not pursued that here. (As the CRCs originate in 
the atomic cores, there is nothing to be gained from using 
Wannier interpolation to calculate them.) Alternatively, 
the Wannier matrix elements could be generated start- 



ing from an all-electron first-principles calculation,"^ in 
which case the present formulation should yield the full 
first-principles orbital magnetization. 

Finally, we mention an issue in all DFT-based studies 
of orbital-ciirrent effects, namely that the accuracy of the 
ordinary exchange-correlation functionals (LSDA, GGA, 
GGA-hU, etc.) has not been well tested in this context. 
A variety of interesting ideas have been proposed for im- 
provcid functionals, ■'^^^^ but c;xploring these would take 
us outside the scope of the present work. 



VI. RESULTS 

In this section, the Wannier interpolation method is 
used to calculate the orbital magnetization of the ferro- 
magnetic transition metals bcc Fe, hep Co, and fee Ni. 
We begin by carrying out convergence tests with respect 
to BZ sampling. Converged values are then tabulated for 
the three materials, and compared with measurements 
and previous calculations. Finally, we investigate how 
M is distributed in k-space in the case of bee Fe. 



A. Convergence studies 

Recall that two separate BZ grids are employed at dif- 
ferent stages of the calculation (the ab initio grid used to 
evaluate the Wannier matrix elements, and the interpola- 
tion grid used to carry out the BZ integrals in the orbital 
magnetization expression) , and both must be checked for 
convergence. 

Tabic I shows the calculated orbital magnetization as a 
function of the number N x N x N of points on a uniform 
interpolation grid in the BZ, for a fixed 10 x 10 x 10 
ab initio grid. For A'^ = 20 the orbital magnetization 
per atom is already reasonably well- converged (to within 
0.002 ^b) in the case of Fe and Co, while Ni requires N = 
50 to reach a similar level of convergence. Setting A'' = 
100 allows to converge M to better than 0.0002 /zs/atom 
across the board. By comparison, the calculation of the 
AHC converges much more slowly.^^ With A^ = 25, for 
example, the AHC of bcc Fe is axy = 554 S/cm, about 
73% of the converged value of 756 S/cm, which demands 
a nominal mesh of the order of A^ = 200, adaptively- 
refined around the strongest Berry curvature spikes. 

Next we look at the convergence properties with re- 
spect to the ab initio mesh, keeping the interpolation 
mesh fixed at = 125 (Table II). The situation is now 
reversed, with the orbital magnetization converging rela- 
tively slowly compared to the exponentially fast conver- 
gence reported in Ref. 25 for the AHC. The term Ml'~^ 
actually converges very rapidly, like the AHC, but the 
convergence rate of is held back by the larger term 

In order to shed light on this behavior, we show in Ta- 
ble III the breakdown of M^^ and M^'^, calculated from 
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TABLE I. Convergence of the orbital magnetization of bcc Fe, 
hep Co, and fee Ni (in units of /is /atom) with respect to the 
interpolation mesh in the Brillouin zone, for a fixed 10 x 10 x 10 

ab initio mesh. For each material the magnetization is along 

the easy axis (see Table IV below). 



Interpolation mesh 


bcc Fe 


hep Co 


fee Ni 


10x10x10 


0.0769 


0.0900 


0.0461 


15x15x15 


0.0797 


0.0839 


0.0394 


20x20x20 


0.0731 


0.0830 


0.0455 


25x25x25 


0.0748 


0.0827 


0.0535 


50x50x50 


0.0749 


0.0840 


0.0462 


75x75x75 


0.0760 


0.0841 


0.0472 


100x100x100 


0.0761 


0.0838 


0.0466 


125x125x125 


0.0760 


0.0839 


0.0468 


150x150x150 


0.0760 


0.0840 


0.0468 


TABLE II. Converg 


;encc of the orbital 


magnetization Mz (in 


units of /is/atom) of bcc Fe with respect to the ab initio 


mesh, for a fixed 125x125x125 interpolation mesh. The two 


gauge-invariant contributions to = 




are also 


shown. 








Ab initio mobli 


Total M,_ 




All"' 


4x4x4 


0.0855 


0.1050 


-0.0195 


6x6x6 


0.0784 


0.0970 


-0.0186 


8x8x8 


0.0765 


0.0948 


-0.0183 


10x10x10 


0.0760 


0.0943 


-0.0183 


12x12x12 


0.0760 


0.0943 


-0.0183 



Eqs. (71) and (72), into the three types of terms intro- 
duced in Eq. (60). The J2 terms give by far the largest 
contribution to M]^(k), similar to what was found pre- 
viously for the Berry curvature. This is not, however, 
the case for M_^^(k), where the JO and Jl terms make 
comparable contributions,^^ and these terms are the ones 
limiting the convergence rate. The reason is that they 
depend on matrix elements involving the position oper- 
ator, Eqs. (82)-(84). In our implementation such ma- 
trix elements are evaluated on the ab initio grid using 
the finite-differences expression (85), and this introduces 
a discretization error which decreases slowly with the 
grid spacing b. Instead, the J2 terms depend exclusively 
on the Hamiltonian matrix elements (81), whose conver- 
gence rate is only limited by the decay properties of the 
WFs in real space (it is therefore exponentially fast). It 
should be possible to achieve an exponential convergence 
for the matrix elements (82) (84) by evaluating them di- 
rectly on a real-space grid, but we have not explored that 
possibility in our calculations. 

Based on the results of the convergence tests presented 
here, we ultimately chose to work with a 10 x 10 x 10 ab 
initio grid and a 100 x 100 x 100 interpolation grid for all 
the calculations presented in the following section. This 
choice of parameters ensures that the values reported for 



TABLE III. Decomposition of the 10x10x10 row of Table II 
into the three types of terms appearing in Eqs. (71) and (72), 
classified according to the number of occurrences of the ma- 
trices [Eq. (76)]. 





JO 


Jl 


J2 




0.0397 


0.0250 


0.0296 




-0.0002 


0.0023 


-0.0204 



the orbital magnetization are converged with respect to 
k-point sampling to within 0.0002 jUs/atom. 



B. Orbital magnetization of Fe, Co, and Ni 

For each of the three materials, two separate sets of 
calculations were carried out, one with the spin magne- 
tization pointing along the easy axis and another with 
the magnetization constrained to point along a different 
high-symmetry direction. In each case the calculated or- 
bital magnetization was formd to be parallel to the spin 
magnetization, as expected from symmetry. 

The numerical results are summarized in Table IV, 
where they are compared with measurements and previ- 
ous calculations. In view of the uncertainties in the accu- 
racy of ordinary GGA functionals for describing orbital- 
current effects, as mentioned at the end of Sec. V, the 
overall agreement with experiment is quite reasonable. 
It can be seen that calculations based on Eq. (1) (both 
ours and those of Ref. 14) give the ordering fee Ni < 
bcc Fe < hep Co for the orbital magnetization per atom, 
in agreement with experiment.^ Instead, the approximate 
muffin-tin scheme switches the first two, because of a 
large contribution in bcc Fe coming from the intersti- 
tial regions between the muffin-tin spheres. The calcu- 
lated anisotropy (orientation dependence) of M is very 
small, and agrees reasonably well with the one obtained 
in Ref. 14. 

While they agree in the general trends, some differ- 
ence can be seen between the values of M obtained 
from Eq. (1) of this work and in Ref. 14. Those dif- 
ferences can probably be attributed to a combination of 
several technical factors, including differences in pseu- 
dopotentials, k-point sampling, and our neglect of the 
core-reconstruction corrections. 

Regarding the two gauge-invariant contributions to M 
in Eq. (3), we find that they have opposite signs in bcc 
Fe and hep Co, and the same sign in fee Ni. In bcc Fe 
and hep Co jM'^^l is larger than |M^'-'| by a factor of 
about 5, while in fee Ni that factor is more than 15. 



C. Distribution of orbital magnetization in k-space 

In order to understand in more detail the results of 
Sec. VIA for the convergence of the orbital magnetiza- 
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TABLE IV. Orbital magnetization (in units of /is /atom) for 
bulk Fe, Co, and Ni. Experimental results for M along the 
easy axis were obtained in Ref. 2 by combining saturation- 
magnetization and gyromagnetic measurements. Results from 
the Wannier interpolation of Eq. (1) ("This work"; M^'^ con- 
tributions in parentheses) are given and compared with the 
results obtained in Ref. 14 (see footnote 44) by evaluating 
Eq. (1) without Wannier interpolation ("Modern theory") 
and by integrating currents inside muffin-tin spheres ( "Muffin- 
tin"). All calculations were done using the PBE functional. 







Modern 


theory 


MufRn-tin 


Axis 


Expt. 


Tills work 


K(-f 


ii Ref 14 



bcc Fe [001]^ 

bcc Fe [111] 

hep Co [001] =• 

hep Co [100] 

fee Ni [111]=' 

fee Ni [001] 



0.081 



0.133 



0.053 



0.0761 (0.0943) 
0.0759 (0.0944) 
0.0838 (0.1027) 
0.0829 (0.0999) 
0.0467 (0.0443) 
0.0469 (0.0440) 



0.0658 
0.0660 
0.0957 
0.0867 

0.0519 
0.0556 



0.0433 
0.0444 
0.0868 
0.0799 
0.0511 
0.0409 



*Experimental easy axis. 



tion, let us look at its distribution in k-space, and com- 
pare with the AHC. For the orbital magnetization we 
sum the integrands in Eqs. (4) and (5), 

M,(k) = -2 Im[G,^(k) + H,y{k) - 2£f F,^(k)] , (97) 

and for the AHC wc take the integrand in Eq. (10), i.e., 
the Berry curvature 



n,{k) = -2ImF^y(k). 



(98) 



Wc will examine bcc Fe with the magnetization along 
the easy axis [001], and accordingly we have picked the 
^-components of the axial vectors M(k) and r2(k). 

The; two qiiantities arc plotted in Fig. 2 along the high- 
symmetry lines F-H-P, together with the energy bands 
close to the Fermi level. The Berry curvature (middle 
panel) is notorious for its very sharp peaks, which occur 
when two bands, one occupied the other empty, almost 
touch. ^^'^^ It can be seen in the lower panel that Mz{]s.) 
displays similar sharp features around the same locations, 
but not nearly as pronounced. The reason is that while 
lmGxy(k) and lmHj;y{'k) are individually as spiky as 
fi^(k), a large degree of cancellation occurs when the 
three quantities are assembled in Eq. (97). This explains 
why the calculation of the orbital magnetization is less 
demanding in terms of BZ sampling than the AHC. 

In Fig. 3 we break down the (k) curve of Fig. 2 into 
various parts. The upper panel shows the contributions 
from Eq. (71), i.e., 

Mf (k) = -2Im[if,j,(k) - £FF,^(k)] , (99) 

the middle panel those from Eq. (72), i.e., 

M^C^k) ^ -21m[G,y(k) - SFF,y{k)] , (100) 




P 

(1/2, V2, V2) 



FIG. 2. Band structure. Berry curvature [Eq. (98)], and in- 
tegrand of the orbital magnetization [Eq. (97)] calculated by 
Wannier interpolation along the path F-H-P in the Brillouin 
zone. Atomic units (a.u.) are used in the middle and lower 
panels. 



and the lower panel their sum M^(k). Each panel con- 
tains three curves, labeled JO, Jl, and J2 according to 
the notation of Eq. (60) and Table 111. The J2 curves 
arc the most spiky, giving rise to the sharpest features 
in Mz{is.)- This is because the matrices appear twice 
in those terms, making them very sensitive to small en- 
ergy denominators in Eq. (77). The main features we en- 
countered in Table 111 for the integrated quantities can 
already be seen in this figure. In particular, the predom- 
inance of the J2 terms in Ml'^ (k) (note the logarithmic 
scale on the upper panel of Fig. 3), compared to a much 
more even distribution of M^'-^(k) among the three types 
of terms (middle panel) . 



VII. CONCLUSIONS 

We have presented a first-principles scheme, based on 
partially occupied Wannier functions, to efficiently cal- 
culate the orbital magnetization of metals using the for- 
mally correct definition for periodic crystals, Eq. (1). 
The localization of the WFs in real space is exploited to 
carry out the necessary Brillouin-zone integrals by Wan- 
nier interpolation, starting from the real-space matrix el- 
ements of a small set of operators [Eqs. (74) and (78-80)]. 
The same type of strategy has previously been used to 
evaluate other quantities, e.g., the anomalous Hall con- 
ductivity^^ and the electron-phonon coupling matrix ele- 
ments,'*'' which are notoriously difficult to converge with 
respect to k-point sampling. 

As a first application, we used the method to calcu- 
late the spontaneous orbital magnetization of the bulk 
ferromagnetic transition metals. Compared to the AHC 
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N o 



-4 

O 

-6 




JO terms 
J1 terms 
J2 terms 




r 

(0, 0, 0) 



H 

(1,0, 0) 



p 

(1/2, 1/2, 1/2) 



FIG. 3. Decomposition of the quantity M^(k) = Ml''^{k) + 
Ml'-' (k) plotted in the lower panel of Fig. 2 into contributions 
from the three types of terms defined in Table III. 



in these systems, we find that the orbital magnetization, 
while displaying similar spiky features in k-spaee around 
the Fermi surface, is somewhat less demanding. Never- 
theless, well-converged results still require a fairly dense 
sampling of the BZ, making it advantageous to use an 
accurate interpolation scheme instead of a direct first- 
principles calculation for every integration point. 

Besides being computationally very efficient, the Wan- 
nier interpolation approach has the appealing feature 
that the evaluation of Eq. (1) is done outside the first- 
principles code in a post-processing step. The algorithm 
is in fact completely independent of such details as the 
basis set used in the first-principles calculation. As a 
result, the calculation of the orbital magnetization only 
needs to be coded once in the Wannier package, and the 
same code can be reused by interfacing it with any de- 
sired k-space electronic structure code. 

We envision several possible applications for the 
method developed in this work. One example is the 
study of enhanced orbital moments in low-dimensional 
systems, such as magnetic nanowires deposited on metal 
surfaces. ^'^^ The method is not restricted to the spon- 
taneous orbital magnetization in ferromagnets; changes 
in magnetization induced by perturbations that preserve 
the lattice periodicity can also be treated within the same 
framework. One such application is the determination of 
the NMR shielding tensors using the so-called "converse" 
approach. ^^'^^ At present the converse approach, using 



Eq. (1) for the induced orbital magnetization, has only 
been applied to molecules and bulk insulators. However, 
by combining it with the present Wannier-based formula- 
tion, it could provide a practical route for the evaluation 
of the shielding tensors in metals, a problem which is 
known to demand a very dense sampling of the Brillouin 



zone 
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Appendix A: Derivation of Eq. (64) 

With the help of Eq. (23), the left-hand-side of Eq. (64) 
can be written as 



tr [/H (U^Ao^U) /H (U^BffU)] 
Next we use the identities 

U^A^U = A« - J« 

and 

U^BpU 



(Al) 
(A2) 
(A3) 



(these can be obtained by inserting Eq. (15) into the 
definitions (28) and (33)), to recast Eq. (Al) as 



Now we note that 



(A5) 



since by construction the occupied Hamiltonian-gauge 
states are eigenstates of the Hamiltonian (see Sec. II B). 
Using Eq. (A5) in the expression (A4), factoring out £m, 
and invoking Eq. (A2) once more, yields 

tr [/« {U^AJJ) /"H« {U^ApU)] . (A6) 

Inserting \ = U'^U between and H^, the right- hand- 
side of Eq. (64) is finally obtained. 
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